/*---------------------------------------------------------------------------*\
License
    This file is not part of OpenFOAM but a modified version of original 
    routine distributed with OpenFOAM.

    You can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation; either version 2 of the License, or (at your
    option) any later version.

    OpenFOAM is a trademark of OpenCFD.

    See GNU General Public License at <http://www.gnu.org/licenses/>

Application
    wallShearStressRho

Description
    Calculates and writes the wall shear stress, for the specified times.

    ATS4i modified original code in order to calculate the wall heat
    flux for RAS compressible flow since the original version was only
    incompressible.

    This version is not endorsed by orginal authors. Use at your own risk.

Authors
    Original code by OpenCFD
    Modified version by ATS4i - Dec/2010

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
//#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
//#include "RASModel.H"

#include "basicPsiThermo.H"
#include "compressible/RAS/RASModel/RASModel.H"
//#include "mutWallFunction/mutWallFunctionFvPatchScalarField.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    timeSelector::addOptions();
    #include "setRootCase.H"
    #include "createTime.H"
    instantList timeDirs = timeSelector::select0(runTime, args);
    #include "createMesh.H"

    forAll(timeDirs, timeI)
    {
        runTime.setTime(timeDirs[timeI], timeI);
        Info<< "Time = " << runTime.timeName() << endl;
        mesh.readUpdate();

//        #include "createFields.H"

    Info<< "Reading field U\n" << endl;
    volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

   IOobject rhoHeader
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::NO_WRITE
    );

    if (!rhoHeader.headerOk())
    {
        Info<< "    no rho field" << endl;
    }
    else
    {

    Info << "Reading field rho\n" << endl;
    volScalarField rho(rhoHeader, mesh);

    #include "compressibleCreatePhi.H"

    autoPtr<basicPsiThermo> pThermo
    (
        basicPsiThermo::New(mesh)
    );

    basicPsiThermo& thermo = pThermo();

    autoPtr<compressible::RASModel> RASModel
    (
        compressible::RASModel::New(rho, U, phi, thermo)
    );


        volSymmTensorField Reff(RASModel->devRhoReff());

        volVectorField wallShearStress
        (
            IOobject
            (
                "wallShearStress",
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            mesh,
            dimensionedVector
            (
                "wallShearStress",
                Reff.dimensions(),
                vector::zero
            )
        );

        forAll(wallShearStress.boundaryField(), patchi)
        {
            wallShearStress.boundaryField()[patchi] =
            (
                -mesh.Sf().boundaryField()[patchi]
                /mesh.magSf().boundaryField()[patchi]
            ) & Reff.boundaryField()[patchi];
        }

        wallShearStress.write();
	}
    }

    Info<< "End" << endl;

    return 0;
}


// ************************************************************************* //
